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We study the generation of entanglement between two distant qubits mediated by the surface 
plasmons of a metallic waveguide. We show that a V-shaped channel milled in a flat metallic surface 
is much more efficient for this purpose than a metallic cylinder. The role of the misalignments of 
the dipole moments of the qubits, an aspect of great importance for experimental implementations, 
is also studied. A careful analysis of the quantum-dynamics of the system by means of a master 
equation shows that two-qubit entanglement generation is essentially due to the dissipative part of 
the effective qubit-qubit coupling provided by the surface plasmons. The influence of a coherent 
external pumping, needed to achieve a steady state entanglement, is discussed. Finally, we pay 
attention to the question of how to get information experimentally on the degree of entanglement 
achieved in the system. 

PACS numbers: 03.67.Bg, 42.50.Ex, 42.79.Gn, 73.20.Mf 



I. INTRODUCTION 

In the last years an intense effort has been made to 
control and tailor the coupling between quantum emit- 
ters and the electromagnetic (EM) field. One major force 
driving the interest in this research area lies in Quan- 
tum Information science, which often requires the trans- 
fer of quantum states between matter and light degrees 
of freedom^. Applications such as quantum teleporta- 
tion, quantum cryptography, and other two-qubit opera- 
tions are additionally based on the availability of entan- 
gled two-qubit systems. There have been many works 
analyzing the coupling of qubits provided by the inter- 
change of fermions or bosons^ ^ and, in particular, ad- 
dressing the generation of entanglement due to the cou- 
pling to a common battPM^. Within this context, the 
EM field may constitute the agent needed to prepare a 
system in a targeted entangled state or to couple two 
preexisting entangled systems. In a number of proposed 
schemes qubit-qubit interactions are provided either by 
coupling to a common EM cavity mode^ or, when large 
separations between the qubits are desired, to a guided 
modd^^^^l With independence of the chosen arrange- 
ment, the dominant paradigm in quantum-state engineer- 
ing relies almost exclusively on exploiting the coherent 
dynamics in order to implement t he op erations needed 
for quantum information processing^^l^. The traditional 
view holds that dissipation, being responsible for deco- 
herence, plays only a negative role. However, it has been 
recently realized that the dissipative dynamics associated 
with the coupling of the system to external reservoirs 
can be engineered in such a way that it can drive the 
system to a desir ed sta te encoding the output of a quan- 
tum computatiorP^l^. Implementation of such ideas has 



shown their tremendous potential demonstrating, among 
other results, the generation o f enta ngled states both in 
theor>*^^'^^ and experimentall}*^^^. 




FIG. 1: (Color online) Two qubits separated a horizontal 
distance d are positioned at a vertical distance h from the 
bottom of a channel waveguide milled in a metallic surface. 
The plasmon modes supported by the structure mediate the 
electromagnetic interaction between the qubits. 

Many structures have been proposed to increase 
light-matter interaction, including photonic crystal 
cavities^^ and waveguides^ photonic nanowireJ^, 
and dielectric slot waveguides^. A crucial requirement 
for such devices is the enhancement of the EM field, 
leading to a large Purcell factor, defined as the decay 
rate of the emitter in the presence of the structure nor- 
malized to the decay rate in vacuum. Electric field in- 
tensification is favored by a tighter confinement of the 
EM modes. In connection with this, metallic structures 
are known to support surface plasmon modes propagat- 



ing at the metal interface and displaying strong field 
concentratioiP^. This modal confinement can reach even 
the subwavelength level^^ , a feature extensively exploited 
in plasmonics, e.g.^ for dense waveguide integration^^. 
The interaction with surface plasmons has been also em- 
ployed to control certain properties of quantum emitters, 
including the decay rate^^, angular directionality^, and 
energy t ransfe J^^^. Single plasmon generation^^ and 
detectio d^^ l ^^ l have been experimentally demonstrated, 
and the achievements on plasmon transport switching^ 
and plasmon-assisted qubit-qubit interactioi]^, suggest 
the on-chip implementation of quantum operations in- 
volving qubits in a plasmonic waveguide network. Along 
this line, we have recently explored the generation of en- 
tanglement between two qubits linked by a plasmonic 
waveguide (PW) consisting of a V-shaped channel milled 
in a flat metallic surface and operating in the optical 
regim d^^ ' ^^' (Fig. [l]). In our previous work, we showed 
that the mentioned configuration enables the sponta- 
neous formation of a high degree of entanglement, even 
for qubit-qubit separations larger than the operating 
wavelength. In the present paper a more detailed anal- 
ysis of the two-qubit entanglement generation mediated 
by plasmons is provided, emphasizing its essential rela- 
tionship with the dissipative character of the effective 
two-qubit dynamics. In addition, a more systematic ex- 
position of several aspects of this problem is presented. 
First, we consider two different waveguide geometries, 
cylindrical and channel-shaped, analyzing the impact of 
the waveguide type on the attainable entanglement de- 
gree. The role of dipole moment misalignments is also 
assessed, which is of great importance for experimental 
implementations due to the difficulty in the controlled 
emitter positioning. The influence of the intensity of the 
coherent external pumping, needed to achieve a steady 
state entanglement, is discussed. We also pay attention 
to the question of how to detect experimentally the de- 
gree of entanglement achieved in the system by measur- 
ing cross terms of a second order coherence function. Fi- 
nally, we study the effect of pure dephasing produced by 
non-radiative mechanisms. 

The rest of the paper is organized in five sections: 
Sec. II contains a short description of the setup and 
the main PWs characteristics. In Sec. Ill we summarize 
briefly the formalism of the master equation governing 
the effective two-qubit dynamics. In addition, the mea- 
sures of entanglement and correlation are recalled. Sec- 
tion IV describes how the classical Green's tensor and 
the associated coupling constants entering the master 
equation are computed. Then, the influence of various 
aspects such as waveguide type, emitter position, and 
dipole moment orientation are analyzed. Once these re- 
sults are available, the generation of entanglement with 
or without external laser pumping is discussed in Sec. V, 
and its relation with the dissipative dynamics is high- 
lighted. We also study the relationship between entangle- 
ment and photon-photon correlations and the influence 
of the pumping rate and pure dephasing. Section VI is 
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FIG. 2: (Color online) Transverse cross section of a cylindri- 
cal nanowire (a) and a channel waveguide (b). The color scale 
in (a) and (b) renders the transverse electric field amplitude 
of the supported plasmonic modes, and the arrows show the 
electric field polarization, (c) Dispersion relation for the fun- 
damental mode of the cylindrical nanowire (black circles) and 
channel waveguide (red triangles). The vacuum light line, the 
dispersion relation of a plasmon on a flat silver surface, and 
that of the second mode supported by the cylinder are also 
plotted. 



devoted to the conclusions. 



II. SYSTEM DESCRIPTION AND PLASMONIC 
WAVEGUIDES CHARACTERISTICS 

The system analyzed in this paper consists of two iden- 
tical quantum emitters positioned in closed proximity to 
a metallic waveguide (Fig.jl]), in such a way that their EM 
interaction is dominated by the plasmonic modes sup- 
ported by the quasi one-dimensional structure. The emit- 
ters, which could be atoms, molecules, quantum dots, or 
nitrogen-vacancy centers in diamond, will be modeled as 
two- level systems, with a transition frequency cjq cor- 
responding to an emission wavelength A = 600 nm. A 
point-emitter approach is assumed because it contains 
all the main physics of the problem without involving a 
detailed description of each qubit, which can be cumber- 
some for large molecules or quantum dots^^. In order to 
determine the influence of the PW geometry, we consider 
two different metallic structures: the first is a cylindri- 
cal nanowire and the second a channel waveguide (the 
case depicted in Fig. [T]). These waveguide types have 
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been previously fabricated and successfully demonstrated 
for dense waveguiding^^ and single plasmon generation^^. 
The exact geometry of both structures is detailed in pan- 
els (a) and (b) of Fig. [2] The radius of the cylindrical 
nanowire is R = 35 nm, the depth of the V-shaped groove 
is L = 138 nm, and its angle is ^ = 20°. The consid- 
ered metal is silver, whose electric permittivity at the 
mentioned wavelength is^^ e = -\- isi = — 13 + zO.8. 
The geometric parameters of both structures have been 
chosen so that, at the operating wavelength, only one 
mode is relevant and the propagation length is identi- 
cal for cylinders and channels. The channel waveguide is 
single-moded and the cylinder supports two modes but 
the second one (black dashed line) , being extremely close 
to the light line, is very much extended in the transverse 
cross plane and will not play a relevant role in what fol- 
lows. Since the qubit-qubit interaction will be mediated 
by the plasmonic modes, having identical propagation 
length ensures a meaningful comparison of the results 
obtained with both PWs. The propagation length is 
i = [2/ci]~^ = 1.7 /im, ki being the imaginary part of the 
(complex) modal wave vector, k = k^ ^ ik{. The disper- 
sion relation for both PWs is rendered in Fig. 2|c) and it 
is observed that the curve corresponding to the cylinder 
(black circles) lies to the right of that corresponding to 
the channel (red triangles), implying that the EM field 
of the former is more tightly confined. This is confirmed 
by a comparison of panels (a) and (b), where the trans- 
verse electric field modal profiles and polarizations are 
plotted. For both waveguides the modal size is deep- 
sub wavelength. In spite of the fact that the electric field 
of both structures includes transverse and longitudinal 
components, the former dominate by a factor of about 
10. For this reason it will be later advantageous to orient 
the emitters parallel to the transverse plane. 



III. TWO-QUBIT DYNAMICS, 
ENTANGLEMENT, AND CORRELATION 



In this section the tools required to determine the 
quantum state of two qubits and their entanglement de- 
gree are reviewed. The evolution of the two qubits in 
interaction with the EM field supported by a plasmonic 
waveguide can be represented using a Green's te nsor ap- 
proach to macroscopic quantum electrodynamicfl^^^^'^. 
One important advantage of this method is that all 
magnitudes describing the coupling between the qubits 
and the EM field can be obtained from the classical 
Green's tensor appropriate for the corresponding struc- 
ture. Within this approach, the Hamiltonian for the sys- 
tem in the presence of a dispersive and absorbing material 



is written in the electric dipole approximation as 

J i=l,2 
poo 

— / (ia;[d^E(r^,cj) + h.c.]. 



Here V and f represent the bosonic fields in the medium 
with absorption, which play the role of the fundamental 
variables of the electromagnetic field and the dielectric 
medium. a\ is the i-qubit rising operator, its spatial 
position, uoq is the transition frequency, and "f" stands for 
the adjoint operation. The interaction term includes the 
dipole moment operator = d^(j^+d*(j|, where d^ is the 
dipolar transition matrix element and * denotes complex 
conjugation. In addition, 

E(r,w)=i./^5^ j d3rV£i(r',w)G(r,r',w)f(r',w) 
V Treo J 

(2) 

is the electric field operator. Notice the explicit ap- 
pearance of the Green's tensor G(r, r',a;), which satis- 
fies the classical Maxwell equations for an infinitesimal 
dipole source located at the spatial position r^ Physi- 
cally, the Green's tensor carries the electromagnetic in- 
teraction from the spatial point to r. 

This hamiltonian description is very powerful but, as a 
matter of fact, it contains too much detail for the purpose 
of this paper. The following simpler description, that de- 
rives from the previous one, will be employed here. To 
determine the entanglement of the two qubits induced by 
their EM interaction, we only need an equation governing 
the dynamics of the reduced density matrix p correspond- 
ing to the two-qubit system. Such a representation of the 
dynamics is obtained from Eqs. ([T]) and ([2| by tracing out 
the EM degrees of freedom. The corresponding master 
equation, whose derivation can be found in Refs. I47|48| 
reads as follows: 

dp ir-£r 



alajp-2ajpal), (3) 



where the hamiltonian included in the coherent part of 
the dynamics is 



XI ^^^^ 



a] an. 



(4) 



The interpretation of the various constants appearing in 
Eqs. (|3| and Q is the following. The Lamb shift. Si, is 
due to the qubit EM self-interaction in the presence of 
the PW. At optical frequencies, for qubit-metal distances 
larger than about 10 nm, 5i is very smal l^^^^ * and will be 
neglected in what follows. The level shift induced by 
the dipole-dipole coupling is given by gij, and can be 
evaluated approximately from 



9ij 



heoc' 



rd*ReG(ri,r^-,a;o)dj- 



(5) 
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Finally, the parameters in the dissipative (noncoherent) 
term of Eq. (|3| are given approximately by 

= ^-^d*ImG(ri,r^-,a;o)d^-, (6) 

and represent the decay rates induced by the self (7^^) 
and mutual (7^^) interactions. Expressions ([5| and (|6| 
are obtained by integration of the EM field Green's func- 
tion in the frequency domaiiP^. To reach the result that 
the coherent and incoherent contributions to the coupling 
are proportional to the real and imaginary parts of the 
Green's function, respectively, the Kramers-Kronig rela- 
tion between the real and imaginary parts of the Green's 
function is used^^'^. In deriving the master equation 
a Born-Markov approximation is applied, valid for weak 
qubit-EM field interaction and broadband PWs. Let us 
remark that, as mentioned above, both gij and 7^^ can 
be extracted from the knowledge of and the classical 
Green's tensor in the presence of the PW. The dipole 
moment can be inferred from the measurement of the 
decay rate of one qubit in vacuum, whose Green's tensor 
is well known. Up to this point it has been assumed that 
both dipoles have equal frequencies, but we would like 
to remark that the formalism is a good approximation 
when the frequencies are unequal but sufficiently close to 
each other. In this regard various criteria can be men- 
tioned. According to Refs. 48 and 52 the frequency differ- 
ence should be much smaller than the average frequency, 
whereas Dung and coworkers^'^ state that the frequency 
difference should be smaller than the typical frequency 
scale for which the Green's tensor displays a significant 
variation. We have checked that both criteria are fulfilled 
for dipoles whose emission wavelengths are in the range 
of 600 nm and differ by less than about ten nanometers. 

To solve Eq. (|3| a basis for the vector space corre- 
sponding to the two-qubit system has to be chosen. A 
convenient basis that makes Hs diagonal is formed by 
the following states: |3) = |eie2), |0) = |^i^2), and 
|±) = ^(1^162) ± |ei^2)), where 1^^) (le^)) labels the 
ground (excited) state of the z-qubit. Using this basis 
the evolution of the diagonal elements of Eq. (|3| is given 
by 

P++ W = (7 + 712)P33(^) - (7 + 7l2)p++(^) 
P—{t) = (7-7l2)p33(^) - (7-712)/?— (^) 

Poo{t) = (7 + 7i2)p++(^) + (7-712)/?— (^), 

(7) 

where it has been assumed that the positions and orien- 
tations of the two qubits in their respective planes trans- 
verse to the PW are identical, so that 711 = 722 = 7 and 
7i2 = 721- The diagonal character of in the above 
mentioned basis and the interpretation of Eqs. ^ is de- 
picted in Fig. [3) including the level scheme and the collec- 
tive decay rates induced by the coupling to the EM field. 
Once these decay rates are evaluated in Sec. IV, the gen- 
eration of entanglement will be elucidated with the help 



of this diagram. Notice that the qubit-qubit dissipative 
coupling induces modified collective decay rates (7 + 712) 
and (7 — 712) which, for particular conditions to be de- 
tailed in Sec. V, give rise to subradiant and superradiant 
states. 



|3> 




FIG. 3: (Color online) Scheme of levels for two identical 
qubits located at equivalent positions with respect to the PW 
and with identical orientations (711 — 722 — 7 and 712 — 721). 

Up to now we have assumed that the system evolves 
without the infiuence of any external agent. As a conse- 
quence, the upper levels in Fig. [3] become eventually de- 
populated and the ground level |0), an unentangled state, 
is reached. To prevent this situation, the decays can be 
compensated by externally pumping the two qubits, thus 
maintaining the system in an excited steady state. In 
cavity quantum electrodynamics, the usual situation is 
that of incoherent pumpin^^^I^ due to the practical dif- 
ficulties of coherently exciting qubits which are embedded 
in a cavity. However, our system is geometrically simpler 
and one can produce a coherent pumping by means of a 
laser whose frequency, cjl, is close to resonance with the 
frequency of the qubits^^ The description of this new 
element requires the inclusion of an additional term in 
the hamiltonian of Eq. Q: 

Hl = -\j2ihn,ale"^-*+h.c.]. (8) 

i 

Here the strength and phase of the laser are character- 
ized by the Rabi frequencies fti = d^EL e^^^^^^fi, where 
El and Icl are the amplitude and wave vector of the 
driving laser field, respectively. In the most general case, 
the determination of the density matrix p{t) requires the 
numerical integration of Eq. ^ with appropriate initial 
conditions'^. When the system is pumped, the steady 
state solution can be obtained by setting p = and solv- 
ing the corresponding linear equations. 

In both scenarios (pumped and non-pumped), once 
the density matrix p{t) is known it is possible to com- 
pute various magnitudes of interest, such as those quan- 
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tifying the two-qubit entanglement, or first and second 
order coherence functions, which are directly related to 
measurable properties. Regarding the quantification of 
entanglement, there are several alternatives but all of 
them are related to each other for a bipartite system^^. 
In this paper we make use of the Concurrence^^, which 
ranges from for unentangled states to 1 for maxi- 
mally entangled states, and is defined as follows: C = 
[max{0, V^- V^}], where {Ai, A2, A3, A4} 

are the eigenvalues of the matrix pTp*T in decreasing 
order (the operator T is cFy ^ Gy^ dy being the Pauli 
matrix). Typical measurable magnitudes include two- 
times coherence functions^^ Their calculation is cum- 
bersome but completely standard, since the quantum re- 
gression theorem^^ establishes that any two-times co- 
herence function obeys the same dynamics as that of the 
density matrix i.e., Eq. (|3|. As it will be discussed 
in Sec. V, entanglement is related with coherence func- 
tions at zero delay. These are measurable by means of a 
Hanbury Brown- Twiss- like experiment detecting photon- 
photon correlations in the emission produced by the de- 
excitation of the qubits. One advantage of zero delay 
correlations is that their calculation is simple because it 
does not require the use of the quantum regression theo- 
rem. 



IV. COMPUTATION OF THE GREEN'S 
TENSOR, DECAY RATES, AND 
DIPOLE-DIPOLE SHIFTS 

In this section we compute the Green's tensor corre- 
sponding to the PWs described in Sec. II. This tensor 
encapsulates the infiuence of the inhomogeneous environ- 
ment and is required for the determination of the decay 
rates, 7^^, and dipole-dipole shifts, gij, appearing in the 
master equation. 



approximation, the length / is kept very short in com- 
parison with the emission wavelength (/ = A/330). To 
model infinitely long PWs the spatial domain of inter- 
est is properly terminated with Perfect Matching Layers, 
which absorb the outgoing electromagnetic waves with 
negligible refiection. The size of the simulation domain 
is of the order of 30 A^. A non uniform mesh is employed 
where the typical element sizes are chosen to satisfy the 
following criteria: ~ A/300 in the dipole neighborhood, 
~ A/40 at the waveguide metal interfaces, ~ A/ 12 at 
the planar metal interface surrounding the channel, and 
~ A/4 in vacuum away from the source. 




Vertical distance, h(nm) 



FIG. 4: (Color online) Purcell factor (7/70) versus vertical 
height h of the emitter along the lines displayed in the insets. 
Cylinder (black circles) and channel (red triangles). 



A. Purcell factor 

For very symmetric structures such as metallic 
plane^^ or cylinders^^ analytic expressions for the 
Green's tensor are available, but for the less symmetric 
case of a channel PW numerical simulations are neces- 
sary. Using the relationship!^ 

E(r)=wVoG(r,r')d, (9) 

the Green's tensor can be inferred if the electric field E(r) 
in position r radiated by a classical oscillating electric 
dipole d at the source position is known. We com- 
pute the EM field excited by the dipole source with the 
finite element method (FEM)f^^^ using commercial soft- 
ware (COMSOL). The point dipole is modeled as a linear 
harmonic current of length /, intensity /q, and orienta- 
tion given by the unit vector n. The associated dipole 
moment is^ d = {iIol/uj)n and, to satisfy the dipole 



Following the explained procedure we now evaluate 
Eq. ([6| to compute the total decay rate, 7 = 711, of 
one qubit in the presence of a PW. This magnitude ap- 
pears in Eq. ([7| setting the time scale of the dynamics. 
The Purcell factor, 7/70, is plotted in Fig.|4]as a function 
of the vertical distance h between the PW and the qubit 
along the vertical lines displayed in the insets (70 denotes 
the decay rate in vacuum). To achieve optimal coupling 
the dipole is aligned parallel to the field polarization, 
i.e., vertically for the cylindrical waveguide and horizon- 
tally for the channel. The Purcell factor is strongly en- 
hanced when the emitter is very close to the metal surface 
{h ^ 0). This effect is more pronounced for the channel, 
due to a higher electric field when the emitter lies at the 
bottom of the groove. The curve corresponding to the 
channel waveguide displays distinct oscillations for large 
h. These are the result of constructive and destructive in- 
terference of the direct field and the field refiected mainly 
at the fiat metallic interface surrounding the channel. 
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is not too high. To be more precise, Eq. (10) is the trans- 
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FIG. 5: (Color online) Beta factor (7P1/7) versus vertical 
height h of the emitter along the lines displayed in the insets. 
Cylinder (black circles) and channel (red triangles). 



B. /3 factor 

The total dipole emission that we have just presented 
can be either radiated to vacuum, non-radiativel y ab- 
sorbed in the metal, or coupled to guided modeP^'^. 
It is thus customary to express the total decay rate as 
the sum of those three contributions, 7 = 7r + 7nr + 7pi- 
The photons absorbed in the metal and most photons 
radiated to vacuum do presumably not contribute to the 
qubit-qubit coupling. It will therefore be interesting to 
compute the decay rate to plasmons, 7pi, and the fraction 
of all emission that is coupled to plasmons, /3 = 7pi/7. 
As will be shown later, these magnitudes play a dominant 
role in the qubit-qubit interaction for appropriate qubit- 
PW vertical distance. In a similar way to the above men- 
tioned total decay rate decomposition, the total Green's 
tensor can be separated as the sum of several terms cor- 
responding to the three emission channels. In order to 
compute 7pi, the plasmon contribution to the Green's 
tensor is required, which is given b^^^^ 



Gpi(r,rO 



i E*(r*) (g)E*(r'*) 



Jk{z- 



(10) 



The occurrence of the exponential factor e*^^^"^ ^ mirrors 
the quasi one-dimensional character of the PW-mediated 
interaction. The lateral extension of the plasmon is taken 
into acount by E*(r*) and H*(r*), which are the trans- 
verse EM fields corresponding to the mode supported by 
the PW [Figs. [2|a) and (b) display the transverse elec- 
tric field] and are evaluated at the transverse position 
= (x^y). Soo is the (infinite) transverse area, is a 
longitudinal unit vector, and (g) denotes the tensor prod- 



verse part of the Green's tensor, which is the relevant part 
since we will only consider transversely oriented dipole 
moments. The modal fields entering Eq. (10) are ob- 



tained by FEM numeri cal si mulation of the correspond- 
ing eigenvalue problenP^'^. Inserting Eq. (10) in the 
expression for the decay rate (Eq. |6| we obtain 



pi 



e ''•^^ cos[kr{z 



(11) 

which, for = vj and = dj, becomes the plasmonic 
decay rate, 7pi. This expression clarifies that 7pi is largest 
when the emitter is positioned at the field maximum and 
aligned with the field polarization. Once 7 and 7pi have 
been determined, we can plot the P factor as a function 
of the vertical distance h between the PW and the qubit 
(Fig. 5|. The general behavior is similar for both the 
cylindrical and channel PWs. First, the P factor is very 
low for small emitter-PW distance, in sharp contrast to 
what is observed for the Purcell factor in Fig. |4] The 
explanation is that 7nr behaves as h~^, where h is the 
qubit-metal distance^^, being the dominant contribution 
to 7 for ^ and effectively quenching the plasmon 
emission. For intermediate h the plasmonic decay domi- 
nates and 13 attains a maximum. Finally, for large h the 
emitter is outside the reach of the plasmon mode and the 
unbounded radiative modes have a larger weight leading 
to a decrease in f3. Nevertheless, the precise behavior 
of P is not identical for both PWs. Channels display 
a higher maximum than cylinders (0.91 at = 160 nm 
versus 0.62 at = 20 nm, respectively) and, in addi- 
tion, the maximum is broader for channels than for cylin- 
ders {P deviates less than a 10% of the maximum value 
within a /i-range of Ah = 100 nm for channels and of only 
Ah = 30 nm for cylinders). These features make chan- 
nels a more attractive structure to enhance the interac- 
tion mediated by plasmons, in the range of parameters 
explored. 



C. Dipole-dipole shift and decay rates 

For high P factor, a dipole couples mainly to plasmon 
modes and this, in turn, warrants that the qubit-qubit in- 
teraction is predominantly plasmon-assisted. Under this 
condition, Eqs. ^ and (|6| for gij and jij can be eval- 
uated using the plasmonic contribution of the Green's 
tensor, Gpi(r, r'), of Eq. (10) instead of the total one, 
G(r, r'). The resulting approximations for the dipole- 
dipole shift and decay rates are as follows^^: 



gij ^ gij,p\ = l^Pe '^^'^^ sm{krd) 
lij - 7ij,pi = 7/^e"'^/^^cos(/cr(i), 



(12) 
(13) 



uct. The derivation of Eq. (10) assumes that the mode 
propagates towards the right > z') and its absorption 



where it has been assumed that the transverse position of 
both qubits and their orientations are identical. Notice 
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FIG. 6: (Color online) Comparison of the exact coupling pa- 
rameters {^ij, Qij) with their plasmonic contributions (72^, pi, 
gij, pi) J as a function of the qubit-qubit horizontal separation 
normalized to the plasmon modal wavelength, d/Xpi. All pa- 
rameters are normalized to the vacuum decay rate 70. (a) 
Cylindrical and (b) channel waveguide. The position and ori- 
entations of the dipoles are detailed in the main text. 



that plasmonic decay is accounted for in Eqs. (12) and 
(13) by the presence of the exponential factor e^/^^. In 



order to check the validity of this approximation a com- 
parison of the exact parameters (^^j, jij) and the ap- 
proximate ones (^ij,pi, 7ij, pi) is presented in Fig. [g] for 
the cylinder [panel {a)] and the channel [panel (b)J~^All 
parameters are normalized to the vacuum decay rate 70. 
In both cases, the position and orientation of the qubits 
are chosen to maximize /3, i.e., h = 20 nm and vertical 
orientation for the cylinder, and h = 150 nm and hori- 
zontal orientation for the channel. The parameters are 
represented as a function of the qubit-qubit separation, 
(i, normalized to the modal wavelength, Api = 27r//cr (at 
the operating wavelength Api is 417 nm for the cylinder 
and 474 nm for the channel). As expected, the approx- 
imation is good for the cylinder and excellent for the 
channel, in consonance with the corresponding P factors 



(0.6 and 0.9, respectively). For the cylinder, at the cho- 
sen the radiative modes play a small but non-negligible 
role which shows up as a small disagreement between the 
exact and approximate results. For both PWs and very 
small many radiative and guided modes contribute to 
the interaction and the approximation breaks down. A 
different approach to this issue leading to the same re- 
sult can be found in Ref. |M] The coupling parameters 
gij and 7^^ are functions of the separation d which os- 
cillate with a periodicity given by the plasmonic wave- 
length, Api, and decay exponentially due to the ohmic 
absorption of the plasmonic mode. Notice that the max- 
ima of jij and those of gij are shifted a distance Api/4, 
which implies that the noncoherent and coherent terms 
of the master equation have different weights for different 
qubit-qubit separations, a fact that will be important in 
Sec. V. 



D. Dipoles with different orientations or vertical 
positions 

To close the analysis of the coupling parameters, we 
now discuss the case when the dipoles have different 
dipole moment orientations or vertical locations. This 
is very important from the experimental point of view 
since a controlled positioning of the emitters is techni- 
cally challenging^^^^. When the two dipoles are inequiv- 
alent in orientation or position, the mutual decay rates 
are obtained in a similar way than Eq. (13) and can be 
expressed as 



7ij,pi = ^/Wfjj^/Mj 



-d/2£ 



COs(/Cr(i), 



(14) 



which indicates that /3's and 7's of both dipoles should 
be as high as possible to obtain a high 7ij,pi. Since 
the dependence of f3 with the vertical distance has been 
discussed already, we focus now on the case of identi- 
cal transverse positions but different orientations for the 
dipoles. The dependence of P with the angular deviation 
of the dipole with respect to the electric field polarization 
is illustrated in Fig. [t] Panels (a) and (b) correspond to 
the cylinder and the channel, respectively. In both cases 
the emitter position is chosen to maximize (3 {h = 20 nm 
for the cylinder, and h = 150 nm for the channel). The 
dipole moment is parallel to the transverse plane, and the 
definitions of the angular deviation, 0, are sketched in the 
diagrams of the corresponding panels. As a general rule, 
the deviation of the dipole from the electric field direction 
has a detrimental effect, and P becomes null for = 90°. 
Nevertheless, there is a broad angular range where P re- 
mains relatively stable so that it is not critically affected 
by relatively large misalignements. Figure [7| shows that 
P deviates less than a 10% of the maximum value within 
a ^-range of AO = 60° for cylinders and of = 40° 
for channels. The functional dependence of p with is 
not simple because although 7pi oc cos^ [see Eq. (pTj)], 7 
has a more complex dependence. This can be observed 
by comparison of the curves in the insets of Fig. jT] We 
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and orientations are identical and chosen to maximize 
the /3 factor. In this simple but insightful configura- 
tion, gij vanishes and jij attains its maximum value 
[Fig. [6jb)], which means that the two-qubit dynamics is 
purely dissipative. The system is initialized in the (un- 
entangled) state |1) = \e1g2) = ^(|+) + |— ))• In this 
case the evolution is confined to the subspace spanned by 
5 1+) 5 master equation is reduced to 

= -(7 + 7i2)p++W 

P—{t) = -{1 -ll2)p—{t) 

Pm{t) = (7 + 7i2)p++W + (7-7i2)/>— 
= -7P+-W- 

(15) 

There are only a few non-zero entries in p{t) and the 
resulting expression for the concurrence is very simple: 

C{t) = V[p++W-P-(iF+4Im[p+_(t)]2, (16) 

where we see that an imbalance of the populations p++ 

and p results in a non-zero concurrence {pj^-{t) is real 

for the chosen conditions). Solving Eq. (15), the concur- 
rence becomes 



C{t) = e-^' sinh [j(3e- 



-Api/(2£) 



t]. 



(17) 



This concurrence and the relevant populations are plot- 
ted in Fig. [s] as a function of time {C is the black thick 
line, and p are the red dashed and blue dotted 



FIG. 7: (Color online) Beta factor of one emitter as a function 
of the angle, 0, formed by the electric field and the dipole 
moment, (a) Cylinder, and (b) channel. The insets show the 
total (continuous line) and plasmon (dashed line) decay rates 
normalized to the vacuum decay rate, (7/70, 7pi/7o), as a 
function of 0. The positions of the dipoles are detailed in the 
main text. 



conclude this section with a brief summary of its main re- 
sults. We have derived simplified expressions for gij and 
7ij, which depend on /3 and 7, and the analysis has shown 
that channel PWs display higher values of the later pa- 
rameters. Therefore, to achieve a larger qubit-qubit cou- 
pling, we will mainly focus on channel waveguides in the 
discussion of the generation of entanglement in the next 
section. 



ENTANGLEMENT GENERATION 



0.5 

jO.4 

^0.2 
0.1 

0.0 

0.5i 





(a) 




















A. Spontaneous decay of a single excitation 

We first consider two identical qubits in front of a chan- 
nel PW without external pumping. The qubits sepa- 
ration is set as d = Api and their transverse positions 



FIG. 8: (Color online) Concurrence (black thick line) and 

populations y9++ (red dashed line), p (blue dotted line), 

versus time, (a) Ideal PW satisfying — \ and ^ = cxo. (b) 
Realistic channel PW. The time is scaled with the emitter 
lifetime (I/7). 
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lines, respectively). Panel (a) corresponds to the ide- 
alized case where /3 = 1 and the plasmon propagation 
length is ^ = oo. The entanglement grows with time 
monotonically up to a value of C = 0.5. This process 



can be easily understood using Eq. (16) and observing 
the mentioned population imbalance. Since 712 = 7, the 
population p++ decays at an enhanced rate 27, whereas 
/)__ stays constant due to its zero decay rate. Panel (b) 
corresponds to a realistic channel PW with p = 0.9 and 
I = 1.7 /im. In this case the concurrence reaches a max- 
imum value of C = 0.33 for t ^ I/7 and then decays 
exponentially to zero. Again, the entanglement genera- 
tion is a consequence of the populations imbalance. For 
this realistic structure both populations have finite de- 
cay rates and the concurrence eventually vanishes. The 
same setup with a cylindrical waveguide produces quali- 
tatively similar results as in Fig. [sj^b) but, since (3 = 0.6 
in this case, the maximum of the concurrence is lower, 
C — 0.21. In all three cases, |+) and |— ) are examples 
of superradiant and subradiant states, respectively. We 
can now present a qualitative picture of more general en- 
tanglement generation processes by referring to Fig. |3] 
The upper level depopulates along two routes: through 
the state |+), with decay rate 7 712, and through the 
state |— ) with decay rate 7 — 712. It is the difference 
in the decay rates along both routes what results in the 
transient build up of the concurrence. Notice that the 
magnitude and sign of 712 depend on d (Fig. [6| causing 
that the roles of the states |+) and |— ) are exchanged for 
^ = 1+) being subradiant and |— ) superradiant. 



B. Stationary state under external continuous 
pumping 

We have just seen the spontaneous generation of en- 
tanglement but, as explained above, the process is a tran- 
sient phenomenon. To compensate the depopulation of 
the upper levels, the system could be externally pumped 
by me ans o f a laser in resonance with the frequency of the 
qubit^^^^. The concurrence reached in the correspond- 
ing steady state. Coo, is plotted in Fig. [9] (black lines) 
as a function of the qubits separation normalized to the 
modal wavelength, d/Api. Three kinds of coherent driv- 
ing have been considered, differing in the relative phase of 
the laser fields acting on qubit 1 and 2: symmetric pump- 
ing means identical Rabi frequencies, Vt\ = 1^2 [panel (a)], 
antisymmetric pumping means fti = —Q2 [panel (b)], 
and asymmetric pumping corresponds to l^i 7^ 0, 1^2 = 
[panel (c)]. The absolute value of the non-zero Rabi fre- 
quencies is 0.157 for the asymmetric pumping and O.I7 
for the other two situations, i.e., relatively weak. It is 
very important to realize that we consider now arbitrary 
separations between the qubits and this implies that both 
coherent and dissipative dynamics are active, its rela- 
tive weight depending on d (Fig. [6|. With independence 
of the pumping scheme the concurrences Cqo in Fig. |9] 
present an oscillating behavior with the qubits separa- 




0.6^ 

I— ' 

to to 



FIG. 9: (Color online) Steady state concurrence (black line) 
and qubit-qubit correlation (red dashed line) as a function 
of the normalized separation d/Xpi. (a) Symmetric pumping 
(Qi = = O.I7), (b) antisymmetric pumping (Qi — —Q2 — 
O.I7), and (c) asymmetric pumping (Qi — O.I57, = 0). 



tion, and damped due to the plasmon absorption. Im- 
portantly, the concurrence maxima occur for those (i/Api 
where the absolute value of jij is maximum (Fig. l6|. 
This suggests a relationship between entanglement gen- 
eration and dissipative two-qubit dynamics. Let us jus- 
tify the position of the maxima of Coo applying the ideas 
developed for the undriven case. When the pumping is 
symmetric [panel (a)], the laser populates the symmetric 
state |+). This state is subradiant for d = ^Api, |Api, . . . 
leading to a population imbalance and the corresponding 
concurrence. For d = Api, 2Api, . . ., |+) is superradiant 
and the pumping is not able to induce a significant p++ 
population. For antisymmetric pumping [panel (b)] it is 
the state |— ) which is populated. This state is subradi- 
ant for d = Api,2Api,... again leading to a population 
imbalance and entanglement. For d = ^Api, |Api, . . ., the 
situation is reversed. Finally, for asymmetric pumping 
[panel (c)] both |+) and |— ) are populated and the sit- 
uation is a mixture of the previous two. In this case 
maxima are found for d = |Api, Api, |Api, . . ., their con- 
currence being slightly smaller than that found for the 
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symmetric or antisymmetric pumping. 




ib) 



d/X 



■pi 




|3) 



+ 



FIG. 10: Tomography of the absolute value of the elements 
of the steady state density matrix for asymmetric pumping 
(Qi = 0.157, ^2 = 0). (a) d = Api/2, and (b) d = Api. 

To verify that the previous interpretation is correct, we 
plot the tomography of the steady state density matrix in 



Fig. 10 We choose the case of asymmetric pumping and 
two different qubit separations. In panel (a) d = Api/2 
and, besides the population of the ground state, we recog- 
nize the large p++ population of the subradiant state |+) 
driven by the pumping, and the negligible p popula- 
tion of the superradiant state |— ). For d = Api [panel (b)], 
we now observe a large p population of the subradi- 
ant state | — ) driven by the pumping, and a negligible 
population of the superradiant state |+). Let us 
remark that, strictly speaking, Eq. (16) is not correct 



when pumping is included, because now further elements 
of p are non zero. However, the tomography shows that 
these additional elements are very small and Eq. (16) 



should be approximately valid, justifying the argument 
that population imbalance leads to concurrence. Since 
this population imbalance is due to the different decay 
rates of the super- and subradiant states, both of which 
are produced by dissipation, we want to emphasize that 
the entanglement generation is driven by the two-qubit 
dissipative dynamics. At this point, a brief comparison 



with the results that can be achieved with cavities may 
be useful again. In cavity QED there is mainly coher- 
ent coupling between the qubits but no cross decay, and 
a coherently pumped cavity is unable to generate any 
significant concurrence. It would be possibl e to work 
with an incoherent pumping with cross term4^^^ but, 
as mentioned previously, this scheme is experimentally 
more difficult than our proposal. 

Once the tomography of the density matrix is known, 
the calculation of concurrence (or any other equivalent 
entanglement quantifier) is straightforward. However, 
tomographic procedures are experimentally cumbersome 
and, for this reason, it is of interest to establish connec- 
tions between entanglement and other more easily mea- 
surable magnitudes. In our two-qubit system, entangle- 
ment is associated with the probability that the state of 
the system is |+) or |— ). In other words, entanglement 
is related with having a strong correlation between the 
states |1) = \e1g2) and |2) = \g1e2). This must mani- 
fest in the correlation between one photon emitted from 
qubit 1 and another photon emitted from qubit 2. Han- 
bury Brown- Twiss-like experiments are able of measuring 
photon-photon correlations and, in particular, the cross- 
term of the second order coherence function which, for 
zero delay, takes the formP^l^ 



(2) 
9l2 



< crj(T2Cr2Cri > 
< o^lai >< crla2 > 



(18) 



Figure |9] displays together the concurrence Coo (black 
continuous lines) and the second order correlation func- 

tion at zero delay (red dashed lines). In all three 

panels it is observed that when Coo is large, a clear anti- 

bunching {g\2 0) takes place, which is consistent with 

the system predominantly being in a state |+) or |— ). 

(2) 

On the other hand, when Coo 0, gi2 grows and the 
antibunching is reduced, which is again consistent with 
a decreased correlation between |1) and |2). The main 
result to be drawn is the distinct relationship between 

f2) 

Coo and g\2 . Lacking an analytical expression relating 

(2) 

Coo and gi2 , our results clearly support the idea of mea- 
suring cross-terms of the second order coherence, at zero 
delay, as a manifestation of entanglement. 

Up to now we have considered a weak pumping rate. 
For an experimental implementation of our proposal, it 
is important to determine the pumping rate range for 
which the described phenomena may happen. The in- 
fluence of the pumping intensity is analyzed in Fig. [TTJ 
which renders Coo versus ^^1/7. Here asymmetric pump- 
ing is considered and a qubit separation d = Api/2. The 
results are computed for three waveguides: a cylinder 
{f3 = 0.6, black circles), a channel (/3 = 0.9, red trian- 
gles), and an ideal waveguide (/3 = 1 and no absorption, 
black line). Each structure presents an optimum pump- 
ing power to achieve maximum concurrence. In order 
to obtain a non-negligible concurrence, the subradiant 
state has to be populated at a rate faster than its life- 
time, which explains both why concurrence is small at 
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0.01 



FIG. 11: (Color online) Steady state concurrence as a 
function of the driving laser power for asymmetric pumping 
(Qi 7^ 0, = 0) and qubits separation d — Api/2. Ideal 
case (3—1 (black line), cylinder /3 = 0.6 (black circles), and 
channel /3 = 0.9 (red triangles). 



low pumping rates and why the structures with lower j3 
require a higher pumping to reach their optimum entan- 
glement. In addition, we observe that the maximum at- 
tainable concurrence improves for higher /3 factor, which 
again justifies the use of channel instead of cylindrical 
PWs. 

Finally, we pay attention to how the generation of en- 
tanglement is affected by the presence of dephasing. For 
this purpose we have recomputed the dynamics of the 
system including now in the master equation (|3| an ad- 
ditional term representing pure dephasing. This term is 
given by^^ 



>Cdeph [p] 



E 



(19) 



The value of the dephasing rate ^y^ is difficult to esti- 
mate in general because it is very dependent on the par- 
ticular realization of the qubit and it is strongly influ- 
enced by the presence of the metallic part of the system. 
For nitrogen- vacancy centers in diamond under resonant 
pumping, pure dephasing times up to 100 ns have been 
measurecP^^. For the typically considered situation in 
our system, where the Purcell factor is about 10, this cor- 
responds to about one hundredth of the emission rate 
7. In our calculations we will consider larger dephasing 
values, both as a conservative measure and because they 
may be more relevant for other emitter types. Figure 12 
shows the steady state concurrence as a function of the 
qubit-qubit separation d for different values of the pure 
dephasing rate and various pumping conditions. Dephas- 
ing grows from zero in panel (a) to ^y^/^ = 0.4 in panel 
(c). The qualitative behavior is the same in all panels but 
the value of Coo decreases as the dephasing rate grows 
(notice that the vertical scale is not the same in all pan- 
els). Nevertheless, the value of the concurrence maxima 
are non- negligible even in the worst case of panel (c). 
Moreover, this decrease can be partially compensated by 
increasing the intensity of the pumping laser. Therefore, 




FIG. 12: (Color online) Steady state concurrence as a function 
of the normalized separation d/Xpi for different pumping con- 
ditions and dephasing rates, (a) ^y^ = 0.0, (b) ^y^ = 0.2, 
and (c) 7'^/7 = 0.4. In all panels the blue dotted lines cor- 
respond to symmetric pumping (^^i = = 0-l7)5 the red 
dashed lines correspond to antisymmetric pumping (Qi = 
—^2 = O.I7), and the black continuous lines correspond to 
asymmetric pumping (Qi = O.I57, = 0). Notice that the 
vertical scale is not the same in the three panels. 



our results show that pure dephasing reduces qubit-qubit 
entanglement but not as much as to preclude its forma- 
tion by the mediation of the surface modes supported by 
ID plasmonic waveguides. 



VI. CONCLUSIONS 

We have presented a detailed analysis of how plasmonic 
waveguides can be used to achieve a high degree of en- 
tanglement between two distant qubits. A full account of 
the theoretical framework has been also described. Im- 
portantly, the degrees of freedom associated with the sur- 
face plasmons can be traced out, leading to a master 
equation formalism for the two qubits in which the two 
contributions to the effective interaction between them 
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(coherent and dissipative terms) are then obtained by 
means of the classical electromagnetic Green's function. 
We have shown that the main ingredients to obtain a high 
value for the concurrence are a large /3-factor and the one- 
dimensional character of the surface modes supported by 
the plasmonic waveguide. By studying how steady-state 
entanglement can be generated, we have demonstrated 
that the dissipative part of the qubit-qubit interaction 
mediated by plasmons is the main driving force in or- 
der to achieve entanglement. We have also analyzed the 
sensitivity of this plasmon-mediated entanglement to dif- 
ferent parameters, such as the dipole orientations of the 
qubits, the pumping rate, and the inherent presence of 
dephasing mechanisms in the system. In all cases, we 
have found that the dissipation-driven generation of en- 
tanglement is robust enough to be observed experimen- 
tally by using plasmonic waveguides that are currently 
available. Finally, we have proposed a feasible way to 
measure the emergence of entanglement in these struc- 
tures by establishing a direct link between the concur- 
rence and the cross-term second order coherence func- 



tions that can be extracted from the experiments. We 
would like to emphasize that the scheme presented in this 
paper could be also operative for other types of photonic 
waveguides provided that the two main ingredients de- 
scribed above (large /3-factor and quasi- ID character) are 
present. Our results demonstrate that plasmonic waveg- 
uides can be used as a reliable tool-box for studying and 
devising quantum optics phenomena without the neces- 
sity of a cavity. 
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